function squFun = continuousSquareWave(t, freq, duty, rampRate)
    %  The Continuous Square Wave is based on the following function
    % 
    %                                 1
    %           squFun(t) = ----------------------
    %                        exp(-R*sin(w*t)) + 1
    %
    %  This function crosses 0.5 when sin(w*t) == 0. To realize the "Duty Ratio" characteristic, we can replace
    %  the term "exp(-R*sin(w*t))" with "exp( -R*( sin(w*t+phase)-level ) )". Where level is calculated from
    %  duty ratio and phase is approperiately chosen to ensure that (t=0) is a center of rising edge.
    %
    %
    % [Inputs]
    %   freq:     The frequency of square wave [Hz]
    %   duty:     Duty ratio of square wave. (Default: 0.5 )
    %   rampRate: Slope at the center of rising/falling edge, in unit of [1/s]. (Default: 1e3*freq )
    %
    if nargin < 3
        duty = 0.5;
    end
    if nargin < 4
        rampRate = 1e3*freq;
    end


%     p = inputParser;
%     p.addRequired('freq', @(x)isnumeric(x) && isscalar(x) && (x>0) ); 
%     p.addRequired('duty', @(x)isnumeric(x) && isscalar(x) && (x>0) && (x<1) ); 
%     p.addRequired('rampRate', @(x)isnumeric(x) && isscalar(x) && (x>0) ); 
%     p.parse(freq, duty, rampRate);

    omega = 2*pi*freq;
    level = sin(pi*(0.5 - duty));    % level = 0 @ duty = 0.5;  level = 1 @ duty = 0; level = -1 @ duty = 1
    phase = pi*(0.5 - duty);
    R = 4*rampRate/omega/sqrt(1-level*level);  % based on the slope @ squFun = 0.5
    squFun = 1./(  exp( -R*( sin(omega*t + phase) - level ) )  +  1  );
end